Visualización del conjunto de datos

civica <- read.csv("Puntos_de_Recarga.csv", 
                     header = T, sep=",", 
                     dec = ".")


kable(head(civica)) %>%
  kable_paper("hover", full_width = F)

Gráfica de los datos

library(tmap)

shape_civica <- st_read("Puntos_de_Recarga_C%C3%ADvica.shp")

## Reading layer `Puntos_de_Recarga_C%C3%ADvica' from data source 
##   `D:\Documentos\GitHub\Patrones-puntuales\Trabajo_patrones_puntuales\Puntos_de_Recarga_C%C3%ADvica.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1278 features and 11 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -75.65818 ymin: 6.0507 xmax: -75.20606 ymax: 6.549086
## Geodetic CRS:  WGS 84

region <- st_read("Area-Metropolitana.shp")

## Reading layer `Area-Metropolitana' from data source 
##   `D:\Documentos\GitHub\Patrones-puntuales\Trabajo_patrones_puntuales\Area-Metropolitana.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -75.7601 ymin: 5.9606 xmax: -75.2373 ymax: 6.5072
## Geodetic CRS:  WGS 84

tmap_options(check.and.fix = TRUE, basemaps = "OpenStreetMap")
tmap_mode('view')

tm_shape(region)+
  tm_polygons(alpha = 0.5, border.alpha = 0.7) +
  tm_shape(shp = shape_civica)+
  tm_dots(size = 0.01)

Conjunto de datos en formato Shape file

## read in the data
shape_civica <- st_read("Puntos_de_Recarga_C%C3%ADvica.shp")

## Reading layer `Puntos_de_Recarga_C%C3%ADvica' from data source 
##   `D:\Documentos\GitHub\Patrones-puntuales\Trabajo_patrones_puntuales\Puntos_de_Recarga_C%C3%ADvica.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1278 features and 11 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -75.65818 ymin: 6.0507 xmax: -75.20606 ymax: 6.549086
## Geodetic CRS:  WGS 84

Se trabaja con los municipios del Valle de Áburra exceptuando las localizaciones de los municipios de Barbosa, Girardota y Caldas ya que tienen pocas observaciones y quedan a los extremos generando sesgo en las estimaciones:

# Conversión a tipo factor
shape_civica$municipio <- as.factor(shape_civica$municipio)

# Municipios a eliminar
indices <- shape_civica$municipio %in% c("Barbosa", "Girardota", "Caldas")

# Eliminación de los municipios en el conjunto de dato
shape_civica <- shape_civica[!indices, ]

class(shape_civica)

## [1] "sf"         "data.frame"

Conversión de longitud-latitud a UTM (metros)

# Latitud - Longitud
data_long_lat <- as.matrix(data.frame(x=shape_civica$longitud,
                              y=shape_civica$latitud))

# Datos en metros (UTM)
data_UTM <- project(data_long_lat ,"+proj=utm +zone=18N +ellps=WGS84")

Creación del objeto ppp de las localizaciones de los puntos de recarga de la tarjeta cívica

Se muestra inicialmente el patrón con el borde del Valle de Áburra:

# Lectura del shape de la región de interés:
coordinate.units <- c("metre", "metres")


# region del area metropolitana
region <- st_read("Area-Metropolitana.shp")

## Reading layer `Area-Metropolitana' from data source 
##   `D:\Documentos\GitHub\Patrones-puntuales\Trabajo_patrones_puntuales\Area-Metropolitana.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -75.7601 ymin: 5.9606 xmax: -75.2373 ymax: 6.5072
## Geodetic CRS:  WGS 84

region <- st_transform(region, crs = CRS("+proj=utm +zone=18N +ellps=WGS84"))

# Ventana
win <- as(as(region, 'Spatial'), 'owin')
ppp.locations <- ppp(x = data_UTM[,1],
                     y = data_UTM[,2],
                     window = win)
# Gráfico del objeto ppp
unitname(ppp.locations) <- coordinate.units
ppp.locations <- as.ppp(ppp.locations)
plot(ppp.locations)

• Como hay zonas del Valle de Áburra sin observaciones (zonas rurales) se debe de considerar un nuevo borde o área de estudio para realizar los analisis con el fin de que no queden muchos espacios vacios (las localizaciones estan ubicadas principalmente en zonas urbanas)

Nueva región de estudio (contorno) a trabajar

win <- ripras(ppp.locations)
ppp.locations <- ppp(x = data_UTM[,1],
                     y = data_UTM[,2],
                     window = win)
ppp.locations <- as.ppp(ppp.locations)
plot(ppp.locations)

Estimación de la intensidad

La intensidad es la densidad promedio de puntos o en otras palabras, el número esperado de puntos por unidad de área. La intensidad puede ser constante (uniforme) o puede variar de un lugar a otro (no uniforme o inhomogéneo).

Un enfoque para evaluar la intensidad es dividir el área de estudio en cuadrantes y contar el número de puntos que caen en cada cuadrante. Si los puntos tienen una intensidad uniforme y son completamente aleatorios, entonces los recuentos de cuadrantes deben ser números aleatorios de Poisson con media constante. Podemos probar esa hipótesis usando el estadístico de prueba de bondad de ajuste \(χ^2\)

par(mfrow=c(1,2))
qc.loc <- quadratcount(ppp.locations, nx=3, ny=4)
plot(ppp.locations, pch=3, cex=0.6)
plot(qc.loc, add=T, textargs = list(col='red'))
plot(as.im(ppp.locations,dimyx=3), main="Conteo por cuandrantes")

quadrat.test(qc.loc)

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  
## X2 = 651.8, df = 10, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 11 tiles (irregular windows)

Estimación de la intensidad via Kernel (no parámetrica)

La estimación no parámetrica de la intensidad, se estima mediante de la densidad del kernel:

Hay diferentes formas de estimar la intensidad, uno de ellos es una estimación no paramétrica mediante suavizado del kernel. Si observamos \(n\) puntos \(\{x_i\}^n_{i=1}\) entonces la forma de un estimador de suavizado kernel es:

\[\hat{\lambda(x)}= \frac{1}{h^2} \sum_{i=1}^n \frac{k \left(\frac{||x-x_i||}{h} \right)}{q(||x||)}\]

donde \(x_i \in \{x_1, x_2, ..., x_n\}\) es un punto observado, \(h\) es el ancho de banda, \(q(∥x∥)\) es una corrección de borde para compensar las observaciones que faltan debido a efectos de borde, y \(κ(u)\) es una función kernel bivariada y simétrica.

No existe una regla general para seleccionar el ancho de banda h, que rige el nivel de suavizado. Los anchos de banda pequeños dan como resultado mapas más puntiagudos y los anchos de banda grandes dan como resultado un mapa más uniforme. Es razonable utilizar varios valores según el proceso que se esté considerando y elegir un valor que parezca plausible (Bivand et al. 2008).

• Probando estimaciones de la intensidad via Kernel con diferentes anchos de banda:

par(mfrow=c(2,3), mar=c(0,0,1,2))

# Ancho de banda por defecto
plot(density.ppp(ppp.locations), main = paste("Sigma (defecto):", round(attr(density.ppp(ppp.locations), "sigma"),3)))

# Ancho de banda: diggle
plot(density.ppp(ppp.locations, sigma = bw.diggle(ppp.locations)), main = paste("Sigma (diggle):", round(attr(density.ppp(ppp.locations, sigma = bw.diggle(ppp.locations)), "sigma"),3)))

# Ancho de banda: CvL
plot(density.ppp(ppp.locations, sigma = bw.CvL(ppp.locations)), main = paste("Sigma (CvL):", round(attr(density.ppp(ppp.locations, sigma = bw.CvL(ppp.locations)), "sigma"),3)))


# Ancho de banda: frac
plot(density.ppp(ppp.locations, sigma = bw.frac(ppp.locations)), main = paste("Sigma (frac):", round(attr(density.ppp(ppp.locations, sigma = bw.frac(ppp.locations)), "sigma"),3)))


# Ancho de banda: ppl
plot(density.ppp(ppp.locations, sigma = bw.ppl(ppp.locations)), main = paste("Sigma (ppl):", round(attr(density.ppp(ppp.locations, sigma = bw.ppl(ppp.locations)), "sigma"),3)))


# Ancho de banda: scott.iso
# plot(density.ppp(ppp.locations, sigma = bw.scott.iso(ppp.locations)), main = paste("Sigma (scott.iso):", round(attr(density.ppp(ppp.locations, sigma = bw.scott.iso(ppp.locations)), "sigma"),3)))

k1 <- density(ppp.locations, sigma = bw.scott(ppp.locations))
ancho_banda <- round(bw.scott(ppp.locations), 2)
plot(k1, main = paste("Sigma (scott)", "\n Anchos de banda: ", ancho_banda[1], ",", ancho_banda[2]))

• Después de probar diferentes anchos de banda para la estimación de la intensidad via Kernel, se elije como el mejor ancho de banda el obtenido mediante scott y adicionalmente se trabaja con un ancho de banda escogido manualmente:

scott

intensidad_scott <- density.ppp(ppp.locations, sigma = bw.scott(ppp.locations))


# Estimaciones de la intensidad mediante scott
par(mfrow=c(1,2))
plot(intensidad_scott, main = "Gráfico de probabilidad")
plot(intensidad_scott, main = "Gráfico de probabilidad")
contour(intensidad_scott, add = TRUE)

Manual

intensidad_visual <- density(ppp.locations, sigma = 700)

# Estimaciones de la intensidad visualmente
par(mfrow=c(1,2))
plot(intensidad_visual, main = "Gráfico de probabilidad")
plot(intensidad_visual, main = "Gráfico de probabilidad")
contour(intensidad_visual, add = TRUE)

• Estimaciones de la intensisad en un gráfico 3D:

scott

persp(intensidad_scott)

Manual

persp(intensidad_visual)

Ajuste de modelos para la estimación de la intensidad (parámetrica)

Modelo base:

Un modelo de Poisson homogéneo

Ejemplo: \(λ_θ((x,y))=exp(θ_0)\) se ajusta:

base.loc <- ppm(ppp.locations, ~1)
coef(base.loc)

## log(lambda) 
##   -12.38517

El modelo log-lineal:

Un modelo de Poisson inhomogéneo con una intensidad logarítmica lineal en las coordenadas cartesianas.

Ejemplo: \(λ_θ((x,y))=exp(θ_0+θ_1x+θ_2y)\) se ajusta:

loglin.loc <- ppm(ppp.locations, ~x+y)
coef(loglin.loc)

##   (Intercept)             x             y 
## -1.110952e+01  6.538523e-05 -4.304356e-05

El modelo log-lineal en x:

Un modelo de Poisson inhomogéneo con una intensidad logarítmica lineal en x.

Ejemplo: \(λ_θ((x,y))=exp(θ_0+θ_1x)\) se ajusta:

loglin.loc.x <- ppm(ppp.locations, ~x)
coef(loglin.loc.x)

##   (Intercept)             x 
## -2.745340e+01  3.459933e-05

El modelo log-lineal en y:

Un modelo de Poisson inhomogéneo con una intensidad logarítmica lineal en y.

Ejemplo: \(λ_θ((x,y))=exp(θ_0+θ_1y)\) se ajusta:

loglin.loc.y <- ppm(ppp.locations, ~y)
coef(loglin.loc.y)

##   (Intercept)             y 
##  2.922814e+00 -2.214429e-05

El log-cuadrático con término de interacción

Un modelo de Poisson inhomogéneo con intensidad que es log-cuadrático y tiene un término de interacción en x e y.

\(λ_θ((x,y))=exp(θ_0+θ_1x+θ_2y+θ_3x^2+θ_4y^2+θ_5x⋅y)\) se ajusta mediante

logquad.loc <- ppm(ppp.locations, ~polynom(x, y, 2))
coef(logquad.loc)

##   (Intercept)             x             y        I(x^2)      I(x * y) 
## -6.898678e+03  3.328491e-02 -1.059952e-03 -8.010213e-08  5.291398e-08 
##        I(y^2) 
## -1.593176e-08

El modelo log-cúbico con término de interacción

Un modelo de Poisson inhomogéneo con una intensidad que es log-cubico y tiene un término de interacción en x e y.

Ejemplo: \(λ_θ((x,y))=exp(θ_0+θ1x+θ_2y+θ_3x^2+θ_4x*y+ +θ_5y^2+θ_6x^3+θ_7x^2y+θ_8xy^2+θ_9y^3)\) se ajusta mediante

logcubic.loc <- ppm(ppp.locations, ~polynom(x, y, 3))
coef(logcubic.loc)

##   (Intercept)             x             y        I(x^2)      I(x * y) 
##  3.295078e+05 -2.185390e+00 -6.285938e-02  6.767948e-06 -2.157113e-06 
##        I(y^2)        I(x^3)    I(x^2 * y)    I(x * y^2)        I(y^3) 
##  7.691461e-07 -8.135718e-12  5.467353e-12 -1.843452e-12  8.260726e-15

par(mfrow=c(2,3),mar=c(0,0,1,2))

plot(predict(base.loc), main='Modelo base')
points(ppp.locations, pch=16, cex=0.2)

plot(predict(loglin.loc), main='log-linear')
points(ppp.locations, pch=16, cex=0.2)

plot(predict(loglin.loc.x), main='log-linear en x')
points(ppp.locations, pch=16, cex=0.2)

plot(predict(loglin.loc.y), main='log-linear en y')
points(ppp.locations, pch=16, cex=0.2)

plot(predict(logquad.loc), main='log-cuadrático con interacción')
points(ppp.locations, pch=16, cex=0.2)

plot(predict(logcubic.loc), main='log-cúbico con interacción')
points(ppp.locations, pch=16, cex=0.2)

par(mfrow=c(2,3),mar=c(1,1,2,2))
diagnose.ppm(base.loc, which = "smooth", main = "modelo base")
diagnose.ppm(loglin.loc, which = "smooth", main = 'log-linear')
diagnose.ppm(loglin.loc.x, which = "smooth", main = 'log-linear en x')
diagnose.ppm(loglin.loc.y, which = "smooth", main = 'log-linear en y')
diagnose.ppm(logquad.loc, which = "smooth", main = 'log-cuadrático')
diagnose.ppm(logcubic.loc, which = "smooth", main = 'log-cúbico')

Estimación de la función k asumiendo inhomogeneidad (Kinhom) para el patrón de puntos de recarga de la civica:

plot(Kinhom(ppp.locations), main='clustered', lambda = intensidad_scott)

Estimación de la función g asumiendo inhomogeneidad (pcfinhom) para el patrón de puntos de recarga de la civica:

plot(pcfinhom(ppp.locations), main='clustered', lambda = intensidad_scott)

Envelopes (Monte Carlo tests)

Kinhom

## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.

• Inicialmente queda por encima de las bandas, por lo tanto es notorio el patrón de cluster o agrupamiento.

pcfinhom

## Generating 99 simulations of CSR  ...
## 1, 2,  [etd 4:27] 3,  [etd 4:08] 4,
##  [etd 4:17] 5,  [etd 4:15] 6,  [etd 4:12] 7,  [etd 4:13] 8,
##  [etd 4:03] 9,  [etd 3:59] 10,  [etd 3:54] 11,  [etd 3:49] 12,
##  [etd 3:44] 13,  [etd 3:41] 14,  [etd 3:37] 15,  [etd 3:33] 16,
##  [etd 3:30] 17,  [etd 3:28] 18,  [etd 3:26] 19,  [etd 3:23] 20,
##  [etd 3:20] 21,  [etd 3:19] 22,  [etd 3:16] 23,  [etd 3:14] 24,
##  [etd 3:10] 25,  [etd 3:06] 26,  [etd 3:05] 27,  [etd 3:01] 28,
##  [etd 2:59] 29,  [etd 2:58] 30,  [etd 2:54] 31,  [etd 2:53] 32,
##  [etd 2:50] 33,  [etd 2:47] 34,  [etd 2:45] 35,  [etd 2:43] 36,
##  [etd 2:39] 37,  [etd 2:37] 38,  [etd 2:34] 39,  [etd 2:31] 40,
##  [etd 2:29] 41,  [etd 2:27] 42,  [etd 2:25] 43,  [etd 2:22] 44,
##  [etd 2:19] 45,  [etd 2:17] 46,  [etd 2:14] 47,  [etd 2:12] 48,
##  [etd 2:10] 49,  [etd 2:07] 50,  [etd 2:05] 51,  [etd 2:02] 52,
##  [etd 2:00] 53,  [etd 1:57] 54,  [etd 1:55] 55,  [etd 1:53] 56,
##  [etd 1:50] 57,  [etd 1:49] 58,  [etd 1:46] 59,  [etd 1:43] 60,
##  [etd 1:41] 61,  [etd 1:38] 62,  [etd 1:37] 63,  [etd 1:34] 64,
##  [etd 1:32] 65,  [etd 1:29] 66,  [etd 1:26] 67,  [etd 1:24] 68,
##  [etd 1:21] 69,  [etd 1:19] 70,  [etd 1:16] 71,  [etd 1:14] 72,
##  [etd 1:11] 73,  [etd 1:08] 74,  [etd 1:05] 75,  [etd 1:03] 76,
##  [etd 1:00] 77,  [etd 57 sec] 78,  [etd 55 sec] 79,  [etd 52 sec] 80,
##  [etd 49 sec] 81,  [etd 47 sec] 82,  [etd 44 sec] 83,  [etd 41 sec] 84,
##  [etd 39 sec] 85,  [etd 36 sec] 86,  [etd 33 sec] 87,  [etd 31 sec] 88,
##  [etd 28 sec] 89,  [etd 26 sec] 90,  [etd 23 sec] 91,  [etd 20 sec] 92,
##  [etd 18 sec] 93,  [etd 15 sec] 94,  [etd 13 sec] 95,  [etd 10 sec] 96,
##  [etd 8 sec] 97,  [etd 5 sec] 98,  [etd 3 sec]  99.
## 
## Done.

• La curva estimada toma valores por encima, especialmente para distancia pequeñas toma valores grandes, quedando por fuera y encima de las bandas teóricas, por lo tanto se concluye que es un patrón tipo cluster - agrupamiento.

Conclusión: Tanto con las estimaciones de las funciones k y g, comparadas con la curva teórica de aleatoriedad completa y con el test de Monte Carlo con las envolventes, se concluye que el patrón de las localizaciones de puntos de recarga de la cívica siguen un patrón tipo cluster - agrupamiento inhomogéneo.

• A partir del conocimiento del tipo de patrón se puede observar diversas zonas con mayor probabilidad de presentar puntos de recargas de la tarjeta cívica, además de otras zonas que aunque pertenecen a áreas urbanas no presentan tanto puntos de recarga por ello tienen un mayor potencial para presentar nuevos puntos de recarga de la cívica.

• Con el análisis del patrón puntual, la estrategia es que para una localización de un punto de recarga de la tarjeta civica hay un radio r (aprox 3000) alrededor con alta probabilidad de presentar otro punto de recarga civica, con esto adicionalmente se le puede recomendar a la empresa metro que para las localizaciones de puntos de recarga donde no existan puntos de recarga a menos de un radio r, es recomendable colocar un nuevo punto de recarga de la tarjeta.